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Proper motions of radiative knots in simulations of stellar jets 
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ABSTRACT 

Aims. Elongated jets from young stellar objects typically present a nodular structure, formed by a chain of bright knots of enhanced 
emission with individual proper motions. Though it is generally accepted that internal shocks play an important role in the formation 
and dynamics of such structures, their precise origin and the mechanisms behind the observed proper motions is still a matter of 
debate. Our goal is to study numerically the origin, dynamics, and emission properties of such knots. 

Methods. Axisymmetric simulations are performed with a shock-capturing code for gas dynamics, allowing for molecular, atomic, 
and ionized hydrogen in non-equilibrium concentrations subject to ionization/recombination processes. Radiative losses in [S II] lines 
are computed, and the resulting synthetic emission maps are compared with observations. 

Results. We show that a pattern of regularly spaced internal oblique shocks, characterized by individual proper motions, is generated 
by the pressure gradient between the propagating jet and the time variable external cocoon. In the case of under-expanded, light jets 
the resulting emission knots are found to move downstream with the jet flow, with increasing velocity and decaying brightness toward 
the leading bow shock. This suggests that the basic properties of the knots observed in stellar jets can be reproduced even without 
invoking ad hoc pulsating conditions at the jet inlet, though an interplay between the two scenarios is certainly possible. 
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1. Introduction 
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i Collimated outflows and jets seem to be ubiquitous features 
in astrophysics and are observed over a wide range of spa- 
, tial scales, from several megaparsecs for extragalactic sources 
■ (AGNs) down to a few parsecs for young stellar objects (YSOs). 
' In both cases a chain of bright knots is typically observed along 
the jet, suggesting that similar physical processes may be at work 
in spite of the enormous difference in scale and, most probably, 
, even in the jet composition itself. As far as YSO jets are con- 
1 cerned, the bright knots observed in emission lines (also known 
as Herbig-Haro, HH, objects) represent spectacular tracers for 
these structures over a wide range of wavelengths, and their 
ubiquity suggests that they are a key component of the jet dy- 
namics. 

Though widely obse rved in detailed images and spectra (e.g. 
iReipurth & BaUvl 1200 lb . the precise nature and origin of the 
knots is, however, still the subject of investigation. In the first 
bright section of an optical jet, up to » 0. 1 pc from the source, the 
number of knots ranges typically between 5 and 20. Their pat- 
tern is not static, but moves at a substantial fraction, from 70% 
or more, of the flow speed, a value increasing with distance from 
the source (Eisloffel & Mundt 1992) . On the oth er hand, the 
knots brightness tends to decay along the jet axis (iMorse et al.l 
ll992l:lRav et al.ll996l:lReipurth & Ballvl200ll) . Several jets have 
been imaged at high re solution (^ 0.1") using the Hubble 



iHartigan et aDl2005h . As an example, HST images of HH 46/47 
exhibit a complex structure i n which the [S II J emission decou- 
ples from the Ha emission jHeathcote et alj[l996.) . Both lines 
define a chain of small knots with spacing of 2 - 3". The Ha 
- [S II] difference image shows that the [S II] emission is more 
extended along the flow than in H a emission maps, where fil- 
aments and wisps are visible, produced either at the front of 
internal working surfaces or in shock induced local changes in 
the flow direction. In this framework the H a emission is pro- 
duced in the hotter and more excited thin layer just behind the 
shock front, while [S II] lines come from a denser and cooler 
layer more distant from the front, which has a larger extension 
jBacciotti & Eisloffel|[T999h . HH 47 also shows a set of leading 
bow shocks with spacing much larger than the intervals between 
neighboring knots. 

The current interpretation for the origin of both the beam 
mini bow shocks and the large leading bow shocks in- 
volves the presence of pulsation in the ejection mechanism 
jRaga & KofmanI 119921: IStone & NormanI 119931: [paUe & Ragal 



Space Telescope (HST) jRav et al.lll 996: Reipurth et al. 1997; 
IHartigan et all 120011: iBallv et al.l 12002- Reipurth et al.. .2002: 
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119951: ISuttneret al.1 1 19971: lRagaetal.1 1 19981 l2002h (the work 
ing surfaces scenario). The physical reason for such pul- 
sation, however, has yet not been clearly identified. It is 
possible that som e kind of self-consistent MHD mechanism 
is at work (e.g. lOuved & Pudritzl 1 1 9971) . A pulsating in- 
flow has also been applied to 3-D simulations which include 
the effects of precession due to the rotation of the nozzle 
(Cerqueira & de Gouveia Dal Pino 2004; Cerqueira et al. 200^, 
motivated by the recent observations of toroid al velocities at 
the base of jets associated with T Tauri stars jBacciotti et al.l 
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l200llCoffev et al.ll2004tlWoitas et aUlIOOllCoffev et aT]l2007h . 
Finally, the effects of magnetic fields in axisymmetric simu- 
lations with periodically varyin g inflow c onditions have been 
also studied (Cerqueira et al. 1 99% lo^ ullivan & Rav 2000; 
Stone & Hardee. ,2000; .Massaglia et al.1 12005; ,de Colle & Ragal 



20061) . 



On the ot her hand, othe r HH jets such as HH 30 jRav et al.l 

11996; Bacc iotli et al.lll999l) do not present bow shock-like fea- 
tures in the forbidden lines, and the knots look more like ax- 
ially symmetric blobs aligned with the jet axis. It has been 
suggested in the past that this appearance c an be generated by 
Kelvin-Helmholtz instabilities. In particular, 'Bodo et al.l (1 19941) 
have investigated the growth of Kelvin-Helmoltz instabilities in 
a Cartesian slab of supersonic, adiabatic flow interacting with 
the external matter at rest. In their paper the authors show that 
as these instabilities grow linearly and eventually saturate, a 
diamond-like pattern of internal oblique shocks (lOS) is formed. 
However, the periodic boundary conditions imposed on the lon- 
gitudinal jet propagation direction prevent all effects due to the 
non-stationary jet propagation and necessarily yield lOS patterns 
which are static with respect to the mean flow. 

In the present study we perform numerical simulations of 
cooling jets which originate from a nozzle, to recover the ef- 
fects of non-stationary jet propagation, but we do not impose 
a periodically variable inflow velocity as it is gener ally as- 
sume d. Encouraged by our preliminary results ( R ubini et al.l 
|2004|) . where intermittent patterns of knots were found in some 
particular cases even without ad hoc assumptions on the inflow 
speed, our intention here is to study, in greater detail and for a 
wider choice of parameters, the formation mechanism, kinemat- 
ics, and emission properties of the lOS which are seen to form 
in the region behind the jet head. These arise due to the pres- 
sure gradient between the jet and the external medium, and may 
be ultimately responsible for the observed knotty emission. 3-D 
effects, such as those due to precession motions at the nozzle 
or to kink-like instabilities induced by magnetic field, are ex- 
pected to be more important far away from the source. For this 
reason, a simple axisymmetric hydrodynamical model can be 
used when investigating the region near the nozzle, where knots 
are seen to form and appear still well collimated with the jet 
beam. Moreover, a cooling model based on a three-species net- 
work (neutral, ionized, and molecular hydrogen) can be safely 
assumed since we are not interested in resolving the small cool- 
ing scales behind the leading bow-shock, where ~ 10-20 
spe cies are typically evolved. For recent results on this subject 
see lRaga et al.l (l2007l) and references therein. 

The steady inflow scenario had been investigated by other 
authors, but it was argued that such lOS should be smoothed 
away i n radiative jets by the c ooling losses (e.g. B londin et alj 
Il990l) . ICerqueira et alj ( Il997h also remarked that in magne- 
tized jets the cooling tends to smooth out pinch-like instabil- 
ities. In general, it is agr eed that radiative losses tend to re- 
duce t he sh ock strength (iDownes & Ravi Il998t iMicono et alj 
Il998l I2OOOI) . so the idea has been put aside in favor of a 
pulsating inflow. However, by using the most recent indica- 
tions concerning the physical parameters of the gas at the jet 
base, obtained from spectral diagnostic s of hig h angular reso- 
lution data (e g. rLavallev-Fouquet et all i2000i: iBacciottill2()02l: 



iHartigan et al.ll2004l) . it is possible to determine the initial con- 
ditions of the numerical simulations in a more realistic way than 
has been done in the past. Here we are able to show that lOS are 
able to survive the cooling losses and to form a regular pattern of 
emitting knots on length scales which mostly depend on the pres- 
sure ratio between jet and interstellar medium. Moreover, such 



a pattern is not static, but the individual knots show a degree of 
proper motion even when steady inflow conditions are imposed. 
This result is basically due to the fact that the gas surrounding 
the jet beam, the cocoon, is a highly dynamic and time- varying 
environment. Thus, the interpretation of the optical knots as be- 
ing due to lOS cannot be ruled out and most probably cooperates 
with the other effects successfully proposed so far. 

The paper is structured as follows. In Sect. 2 we discuss the 
equations, the numerical method, and the general mechanism for 
the formation of lOS and for their proper motion. In Sect. 3 we 
present various results from the simulations, including a com- 
parison with observations. Section 4 is devoted to the final dis- 
cussions. 



2. The physical and numerical model 

2.1. Gas dynamic equations and source terms 

In our simulations we solve the system of hydrodynamical fluid 
equations for a cooling gas formed by molecular (H2), atomic 
(HI), and ionized (HII) hydrogen in non-equilibrium concen- 
trations, plus atomic helium (He) and heavier elements in fixed 
concentrations and free electrons (e). The evolution equations 
for such a system, in conservation form as required by shock- 
capturing numerical schemes, are the following: 



dt 



ipu) = 0, 



—(pu) + V ■ (puu + pi) - 0, 

dt 



— ipE) + V ■ (pHu) = ->Srad - Sioa ' ■Sh. , 

dt 



dt 
dt 



+ V-(A?HII«) = A^ion-A^re 



+ V ■ (A^H, m) = -Ndiss- 



(1) 

(2) 
(3) 
(4) 
(5) 



In the above equations p, u and p represent density, velocity and 
pressure for each volume element, / is the identity tensor, E = 
e + ju^ is internal plus kinetic energy per unit mass, and H = 
E + p/p is the specific enthalpy. These evolution equations are 
completed by the relation for the total number density 



Mr 



(6) 



which is used to derive the gas temperature, where A'hii, A^hi, 
Nhj, Ne, and A^He are the number densities for ionized, atomic 
(neutral), molecular hydrogen, electrons, and helium atoms, re- 
spectively. The source terms on the right hand side of Eq. (O, the 
energy equation, are respectively: >S,ad, which contains all radia- 
tive losses excepting those due to H2; >Sh,, where radiative and 
dissociation losses from H2 are included; <Sion, which takes into 
account hydrogen ionization. These contributions, which require 
computation of the particle densities, will be discussed later. 
Note that the non-equilibrium hypothesis leads to the two ex- 
tra continuity equations, namely Eq. (jUl for protons and Eq. (|5]l 
for molecular hydrogen. In the former, the term Nion - A/'iec is 
the proton population variation, per unit of time and volume, 
due to ionization and recombination, whereas in the latter A/diss 
is the molecular hydrogen density decrease rate due to dissocia- 
tion processes. Reformation processes have not been considered 
here, since densities for molecular hydrogen reformation are too 
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small, and reformation time scales are too large compared with 
dissociation time scales. Solar abundances of heavier elements 
in fixed concentrations have been also included, since they are 
important for energy losses though are too small to affect the 
mass. 

For the computation of the total particle density in Eq. (|6]l, 
which is needed to update the temperature field, the following 
steps are required: 

- the total mass density p and the number densities A^e and Na, 
are evolved directly via the respective continuity equations, 
thus are known quantities at each time-step; 

- let us first introduce the total number density of hydrogen 
atoms as 



A^H =A^Hii+A^Hi+2A^H,, 



(7) 



including all possible states. From the relations pe - Num^, 
PHe = 4A^He'Wp (p = Ph + PHe and mp is the proton mass), and 
assuming helium to be present in a fixed concentration with 
PHe/p = 1 /4, we easily find 

A^H = |(p/'«p), A^He = i^(p/mp); (8) 

- the number density of electrons A^e is derived from 



A^e = A^HII + A^era = A^HII + lO^'^A^H, 



(9) 



where A^em is the concentration of free electrons due to met- 
als, again assumed to be a fixed fraction of A^h; 
- the number density of atomic hydrogen A^h i is finally derived 
from Eq. (|7]i. 

Once we have the temperature field, we are ready to compute 
all source terms in the fluid equations. Let us discuss them in de- 
tail. The source term in Eq. (|4|l for ionized hydrogen is provided 
by recombination and ionization rates, respectively given by 

Mec =A?eA^HII«(n Nion^N.NniPiT). (10) 

The functions a and fi of temperature T (in Kelvin degrees) are 
a(r) = 2.O6xlO""r"'/^0 cm^ s"', (11) 



PiT) = 7.80 X 10""r'^^exp(-13.6eV/feBr) cm^ s" 



(12) 



where in Eq. ( fTTI ) is a decreasing function of the temperature, 
ranging approximately from 4 to 1 in the region of interes t (see 
Table 5.2 in Spitzer ( 1978)), whereas Eq. ( fT2b is taken from lLan ^ 
(Il975h for hydrogen neutral atoms mainly in the ground level, as 
expected for a low density gas, and kg is the Boltzmann constant. 
This yields an energy loss contribution >Sion = (Mon - Mec) ■ 
13.6 eV due to ionization of hydrogen in Eq. ([3]l. 

The source term in the molecular hydrogen equation arises 
from molecular dissociation alone, due to collisions with H2, H I, 
HII and electrons. This yields 

Wdiss = A^hJA^H.^H, + A^HI^HI + A^HII^HII + A^e^e], (B) 

in which dissociation rates K-a^ and Km due to collisions 
with neutrals, H 2 and HI, res p ectiv ely, are computed through 
the model by iLepp & Shulll (Il983l) . whereas those due to 
collisions with charg e d part icles, A'hii and A"e, deri v e from 
iHollenbach & McKed (Il989h and iMac Low & Shuill (Il986h . 
Corrections to the critical density in order to separate high and 
low density regimes, as suggested by Martin et al. ( 1996), have 
also been introduced. The corresponding contribution to the total 
energy source is ^Shj = A^d + Ajiss, where Adiss - A/diss -4.48 eV 
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Table 1. Run parameters. According to our definitions, pjet = 
7/0.836 X 10-2"gcm-^ pjet = H 1.381 x lO^'^ergcm-^ The 
jet velocity is always assumed Vjet = 200 km s"'. Case £ is the 
same as D but with restarting inflow conditions. 



comes from dissocia tions while hraj is m ainly due to collisions 
with H I and H2 (see lLepp & S"huil[T983l) . 

Finally, we consider the coolin g function Sgd- For T > 10^ K 
we use the cooling function as in iDalgarno & McCravl ( |1972|) . 
which considers radiative losses due to de-excitation of H I en- 
ergy levels. Here we assume the gas to be optically thin. For 
T < 10"^ K the cooling function is extended by computing en- 
ergy losses mainly due to collisions between electrons and C II, 
O I, N I, Fe II, O II, S II, or, for still lower temperatu res, between 
HI and CH, OL Sill, Fell dBacciotti et al.lll995h . Moreover, 
we must consider that in star-forming regions, refractory species 
concentrate at the surface of dust grains. The cooling func- 
tion has been corrected accordingly to take into account re- 
duced concentrati ons with respect to standard metal abundances 
(ISofia et al.lll994 . 



2.2. Simulation setup and choice of parameters 

The numerical code used in our simulations is a finite volume 
Godunov-type scheme which solves Eqs. ([T]i-(|5]l, in 2-D, by as- 
suming axisymmetry around the jet axis and adopting cylindri- 
cal coordinates (^, r), where f is the distance from the source 
along the jet and r is the radius. We use second order limited 
reconstruction on characteristic variables and an exact Riemann 
solver at cell interfaces to work out numerical fluxes, which in 
our simulations has turned out to be more robust with respect 
to a Roe solver. The jet originates from a nozzle of a given ra- 
dius, located at the left side of the numerical box and then prop- 
agates in the ^ direction. Convergence tests have been performed 
to find the optimal numerical parameters which reconcile accu- 
racy and efficiency. A stretched grid with 150 points is used in 
the radial direction, which uses 20 points to describe the noz- 
zle region, from r = to r = rjet, whereas along ^ we expand 
the numerical box by adding points (up to 200 at most), while 
maintaining a fixed resolution, as the jet propagates in the un- 
perturbed medi um. Details of t he numerical method and tests 
may be found in lLorussol (Il999l) . 

The physical parameters to be initialized at the nozzle and 
in the external, unperturbed interstellar medium (ISM), are 
density, pressure (and, consequently, temperature), ionization 
and molecular hydrogen fractions (defined as ^ xhii = 
A'hii/(A^hi + A^Hii) and xh^ = Nu.JNu, respectively), plus jet 
radius rjet and inflow velocity Vjet. Most of these parameters 
can be derived from observations, either directly or combined 
with spectral diagnostics at moderate and high angular resolu- 
tion ('Bacciot ti & Eisloffell [19991 iLavallev-Fouquet et alj|2000l: 
iBacc iotti et all 20021) . Typical observed tangential velocities are 
about 100-200 km s"' and we assume Vjet = 200 km (which 
corresponds to a Mach number of the order of 20 for the assumed 
temperature, see below). In general, at 100 - 200 AU from the 
source, the ionization fraction x^ and the electron density Ne in 
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the jet are in the range 0.2 - 0.5 and 10^ - 10^ cm"^, respec- 
tively, while the temperature is in the range 6 000 - 10000 K. 
From this point onward the chain of knots begins to be visible 
and the physical quantities mentioned above oscillate around the 
val ues attained at 200 AU from the source (see, e.g.. Figs. 1 and 
2 in lBacciotti et al.iri999l) . According to these analyses we place 
the inlet at about 200 AU from the star 

The precise parameters for all runs are shown in Table [1] 
Rather than giving the jet density and pressure, we fixed only the 
ISM values and we prescribed the two quantities rj - pjet/piSM 
and n = /Jjet/piSM- Namely, the jet density and pressure are de- 
rived from the given values assuming pisM = 0.5p and pisM - P- 
Normalization is against L - lO'*" cm, p = 1.673 x 10"^" g cm"-' 
(corresponding to 10"* hydrogen atoms per cm-'), p = 1.381 x 
10""^ erg cm"-' (corresponding to a temperature of 100 K for a 
purely atomic hydrogen gas at density p). Here we will always 
assume 11 = 60 77, to preserve the expected high temperature ra- 
tio between jet and ISM, with the condition 11 > 1 . The choice of 
n (or 77) plays a key role, since the pressure ratio affects dramat- 
ically the jet behavior and, in particular, the formation of lOS. 
This, in turn, leads to quite different knot patterns, as we show 
in the next subsection. Another parameter affecting the structure 
of lOS, in particular their dimension and spacing, is obviously 
Tjet. Note that cases Ji, S, and C differ only for the density ra- 
tio 77 (and thus also 11), which is decreased from 10 to 0.1 (thus 
we move from heavy to light jets). Run £) is a light jet case 
where parameters have been further optimized to provide knots 
with a structure similar to that observed, whereas run & retains 
the same parameters as in run 3D, but the ejection is turned off 
for a while in order to reproduce some observations where two 
separate chains of knots are found (case of restarting jet, see 
Sect.[X4ll. 

2.3. Formation mechanism of compression regions in 
under-expanded jets 

We now investigate the formation mechanism of compression 
regions and lOS in under-expanded jets, that is when f jet > 
Asm ^ n > 1. Under these conditions, just after the nozzle the 
gas starts to expand downward until the pressure falls below the 
surrounding local pressure in the cocoon environment. Outward 
from this pressure equilibrium point the flow is deviated by the 
inward radial pressure gradient, forming blobs of compressed 
gas aligned along the axis, which will drive the formation of ra- 
diatively emitting knots. This is the formation mechanism of lOS 
in jets, which is well known in gas dynamics. Here we study how 
the resulting expansion and compression features depend on the 
initial parameters. 

To this aim the code has been adapted to simulate the 2-D ex- 
pansion of a circular surface at constant ^, under-expanded with 
respect to the environment, i.e. to the ISM. The circular surface 
mimics an expanding slice of gas perpendicular to the jet axis 
and co-moving with the jet itself. At time f = it represents the 
jet nozzle, which has been initialized with a typical set of inflow 
parameters. In this qualitative model the effects of the longitu- 
dinal propagation of the jet, including re-circulating flows, co- 
coon temporal variability and Kelvin-Helmoltz instabilities aris- 
ing from the longitudinal velocity shear, have been neglected. 
In the present section, both adiabatic (Al, Bl, CI) and radia- 
tive (A2, B2, C2) simulations are performed, with density ratio 
77 = 10,1,0.1 and pressure ratio 11 = 600, 60, 6, respectively. 
The other parameters are the same as in the runs J^, C listed 
in Table [T] respectively, except for the nozzle size, whose radius 
is ten times smaller, 0.4 x 10'"* cm, also leading to a smaller time 



a) " b) 

o 



C 1 2 




Fig. 1. Case Bl. Time evolution of the pressure field. Times are 
f = yr (a), f = 27 yr (b), t = 36 yr (c), f = 51 yr (d). The pres- 
sure scale is /? = 1.381 X 10"'" erg cm^, while here the length 
scale is 10'"* cm. The lower (white) limit in the color bar is 0.5, 
corresponding to the value assumed for the external ISM pres- 
sure. 




36 years 

Fig. 2. Case B 1 . Velocity vectors at f = 36 yr. A shell of gas is 
moving toward the axis pushed by the inward pressure gradient, 
while the external shells are still expanding outward. 



evolution scale. In the jet simulations presented in the following 
sections a more realistic beam width has been used, to match the 
estimated mass loss rate. 

Case Bl (77 = 1, II = 60, adiabatic) is chosen to illustrate the 
general features of the gas radial evolution as a function of time, 
described in Fig.[T] Plot (a) represents the initial condition, with 
uniform pressure distribution inside the jet beam. At f = 27 yr 
(b) the shell is expanding: a strong pressure wave moves out- 
ward, and an evacuated region forms, inside the beam, driven by 
the internal rarefaction wave. When the inward pressure gradient 
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Fig. 3. Plots of the central pressure (in logarithmic scale) versus 
time for the three values of the pressure ratio n - 600, 60, 6 
chosen for our simulations (cases A, B, C, respectively). Both 
adiabatic (1) and radiative (2) runs are shown. 



becomes strong enough, particles from the intermediate shells 
start to move toward the center (c, t - 36 yr). This reverse flow 
can be seen in the velocity pattern of Fig. |2] Particles moving to- 
ward the axis finally generate a central peak of pressure, which 
will result in a radiatively emitting knot (Fig.[TJl). 

The combined effects of radiative losses and different pres- 
sure ratios are shown in Fig. [3] which shows the pressure taken 
on the axis versus time for all cases. Notice that the pressure 
jump, coinciding with the formation of the central peak, does 
not occur at the same time for all runs. In general, the smaller 
the pressure ratio 11, the smaller the knot formation time scale, 
and in turn the smaller the distance from the injection point. In 
the simulation with 11 = 6, case C, a knot quite close to the 
source is produced (at f = 30 yr, corresponding to a distance 
L w 2 X 10'^ cm if the particles co-moved with the jet at the con- 
stant velocity of 200 km s"'). The peak formation is most appar- 
ent for case B, where we clearly see a pressure jump, whereas in 
run A, with the highest pressure ratio n = 600, no knot is gener- 
ated over the length scale corresponding to 1 00 yr, at the same jet 
velocity. As far as radiative losses are concerned, these in gen- 
eral depend on the matter density and tend to delay the pressure 
jump and lower the peak value, as is most apparent by compar- 
ing cases Bl and B2. In the radiative case A2, in particular, the 
jet matter is so overpressured and dense that radiative losses ini- 
tially result in a quick pressure drop larger than for the adiabatic 
case Al . This radiative loss is also able to keep the pressure value 
at the center below the adiabatic value at later times. All cases, 
however, show that the value of the pressure jump is not sub- 
stantially affected by cooling losses. Instead, on-axis peak and 
minimum pressure values are smaller than for adiabatic simula- 
tions. Radiative losses, in fact, effectively lower the gas adiabatic 
index 7, and allow for a denser and cooler central region. 

As already anticipated, the nozzle radius is another key pa- 
rameter that regulates the knots' length scale. The smaller the 
nozzle radius, the shorter the distance traveled by sound waves 
in the time it takes to go from the axis to the beam surface and 
back to the axis, where the compression region forms. We have 
investigated the suitable initial conditions for a jet simulation. 
To summarize, under-expanded light jets with small jet/ambient 
pressure ratio (but still greater than 1), and small nozzle size 
are expected to produce nodular structures more similar to those 
observed, in terms of spacing of the knots and closeness to 



the source. Moreover, contrary to common belief, the formation 
mechanism described in this section allows for proper motion 
of knots in real jets. Unlike in this simple model, where the ex- 
ternal pressure is constant with time, the beam in a stellar jet 
is in fact embedded in a variable cocoon. At a given point ^ on 
the longitudinal axis the local pressure ratio njc(^) between the 
pressure of the jet and the cocoon embedding the beam (which 
is not equal to the constant value of 11) actually changes in time, 
and the knot has to move to adjust the position to match the 
new value of njc(^), according to Fig. [3] Namely, Iljc(^) is ex- 
pected to grow with time, since the pressure inside the beam is 
steadily fed by the nozzle and stays unchanged, while that in the 
cocoon decays because of both cooling losses and lateral expan- 
sion. Therefore, the knot is expected to move downward follow- 
ing the gas stream, at some fraction of the injection speed, as is 
shown in next section. 

3. Simulation results and emission maps 

In the present section we show the results of the (radiative) sim- 
ulations for all the cases listed in Table [T] now restoring the full 
settings for axisymmetric simulations in (^, r). The results are 
illustrated with 2-D maps of the total density on a jet merid- 
ional section and, in some cases, with derived synthetic emission 
maps calculated for the coUisionally excited lines [S II] AA 6716, 
6731. These lines are commonly observed in YSO jets and thus 
the constructed maps represent a good test of the simulation re- 
sults against the observations. In some cases the integrated emis- 
sion from a slice perpendicular to the jet axis (of normalized 
width 1 cm) will be shown as a function of the distance from the 
source, and hereafter will be labeled as E(^). The [S II] emis- 
sion is calculated from the physical parameters determined in 
the simulation by using a public routine for a 5-level collision- 
ally excited atom (A. Raga, priv. comm.), assuming that all S 
atoms are ionized once and adopting a (constant) relative abun- 
dance S/H = 1.6 X 10"^ (for details see .Bacciotti et al..J995,) . 

3.1. Heavy and density-matched jets (cases Ji, S) 

Case the starting model in the parameter space of Table [T] 
represents a heavy (rj = 10), strongly under-expanded jet with 
pressure ratio 11 - 600. This test case does not produce emitting 
knots. The density field (Fig. IH top panel) does not reveal any 
kind of internal structure. Note that in order to allow for a rep- 
resentation of some details of the jet beam, density and emissiv- 
ity maps in the paper have an aspect ratio of the axes' scales far 
from unity. In the bottom panel we report the integrated emission 
function E{^). The figure only shows random discontinuities that 
arise from rings of dense matter in the external cocoon region, 
rather than from blobs of compressed gas on the axis. Due to 
the high pressure ratio at the nozzle, the gas expands from the 
origin to x 150 L, where recompression occurs due to the termi- 
nation shock corresponding to the Mach disk. Velocity, density 
and pressure are rather smooth and uniform from the source to 
die Mach disk (located at ^ 200 Lalt - 1400 yr) and compar- 
isons with observed images of HH objects cannot be attempted. 
Downstream of the Mach disk a secondary jet forms, acceler- 
ated by the local De Laval nozzle generated by a toroidal ring of 
dense matter that fo rms around the Mach disk triple point (e.g. 
iBlondin et al.lll990l) . Such an effect is freque ntly observed in ax - 
isymmetric HD and MHD simulations (e.g. lClarke et aljri986l) . 
while it does not appear in 3-D simulations. 

Case S, the density-matched run with rf - \ and 11 - 60, 
is still far from showing satisfactory features. Density maps are 
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Fig. 4. Case Ji. Top: density map in logarithmic scale. Bottom: 
[S II] emission E(^), in logarithmic scale, integrated on a cylin- 
drical slice (of thickness 1 cm), vs. the distance from the source 
in units of L. 
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Fig. 5. Case S. Density map in logarithmic scale. 



3.2. Light jets (case C) 

In this case we have 77 = 0.1 and 11 = 6, all other parameters 
being unchanged with respect to cases and S. The mass loss 
rate is of the order of 10"^ Mq yr ~', in agreement wit h the esti- 
mates derived from observations (IBacciotti et al.ll2002l) . Figure|6] 
shows the density map (left panel) and a simulated image in the 
light of [S II] lines (right panel), on the same scale and for the 
same time t - 1500 yr A well defined knot appears at ^ « 90 L 
from the source (labeled with A in the synthetic emission map), 
corresponding to the dark throat observed in the density map. 
The estimated luminosity of this single knot is 10^^ erg s"'. 
This case confirms that lOS and emission knots appear for small 
n values. The simple model that has been used to show how 
blobs of compressed gas form in under-expanded jets, can now 
be checked. Figure [T] contains a sequence of plots showing the 
radial profiles p{r) of the pressure, taken at different distances 
^ along the jet axis. In each plot the solid line refers to time 
t = 900 yr, while the dotted line refers to the same output time 
inPig.ilf = 1500 yr 

We will now discuss the figures in some detail. Panel (a) of 
Fig. |2]refers to a region quite close to the source, ^ = 0.35 L. 
The pressure on the axis corresponds to the inflow value that 
fills the nozzle homogeneously, the external flat profile matches 
the ISM pressure. The two plots show that close to the source 
the pattern does not change in time significantly. Panel (b) gives 
the situation approximately 20 nozzle radii downstream of the 
source. As particles flow down the nozzle, the beam expands 
laterally, following the outgoing pressure wave (see Fig. [TJ) for 
a comparison with the slice model). Comparison between solid 
and dotted line in panel (b) shows that the beam undergoes a lat- 
eral expansion in time. At f = 1500 yr the beam is wider and 
the average cocoon pressure is lower with respect to f = 900 yr 
(the cocoon is defined here as the region bounded by the peak of 
outgoing pressure wave). Figure |7};, to be compared to Fig.[T]:, 
shows the growth of the inward pressure wave that pushes the 
gas toward the axis. This inward flow is plotted in Fig. [8] which 
shows the radial profile of the radial velocity Vr{r), taken at the 
age of f = 900 yr and for ^ - 11 L (quite near the source, solid 
line) and for ^ = 21 L, (closer to knot A, dotted line). The for- 
mer is overall positive (expansion phase), the latter shows that 
internal rings of matter are moving to the axis (compare also 
with Fig. |2]i. In panel (d) the internal compression wave is still 
moving inward, until it reaches the axis (e, see also Fig.[T]i). The 
peak of pressure on the axis in panel (e) corresponds to knot A 
in the young jet (i.e. that at f = 900 yr). The oW jet wave is still 
on the way and reaches the axis in panel (f), approximatively 
40 nozzle radii downstream. This spatial delay confirms that as 
time goes by the knot forms farther and farther away from the 
source, because of the lower external pressure field, according 
to Fig. |3] The spatial gap between the peaks in panels (e) and (f) 
divided by the time interval between what we have labeled as old 
and young jets measures the knot proper motion. In this case the 
resulting velocity is too small when compared to observations 
(Vkno, ^ 9 km s-')- 



similar to those of case as shown in Fig. |5] Although the 
jet is more coUimated than in case y[, no nodular structure is 
visible, since the pressure ratio 11 at the nozzle is still too large 
to produce knots over the jet length. Note also that in this case 
the jet presents a narrower nose cone ahead of the Mach disk. 



3.3. Ligint jets w/fh smaller radius (case D) 

To allow for the formation of knots over the jet length scale and 
to reduce the intra-knot spacing, in case 2D the nozzle radius rjet 
has been reduced to 0. 1 L, while t] has been increased to 0.4, pre- 
serving the usual inflow speed. The increase in density is needed 
to keep a realistic mass injection rate, by (partially) compensat- 
ing the loss of area. A spectacular chain of knots appears, very 
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Log density [g cm^^J at t=1500 yr [SHJ log emission mop [erg s" ' cm"3j at t=1500 yr 
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Fig. 6. Case C. Left panel: density map in logarithmic scale. Right panel: [S II] emission map in logarithmic scale. Note the presence 
of a knot (A). Its estimated luminosity is ~ 10^^ erg s~^ 
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Fig. 7. Case C, Pressure p(r) as a function of the radius, at f = 900 yr (solid line) and t = 1500 yr (dotted line). Panels (a-f) refer 
respectively to positions ^ = 0.35 L, ^ - \ \ L, ^ - 2\ L, ^ - 67 L, ^ - 74 L, ^ = 91 L. Pressure and radius scaling changes from 
picture to picture, to allow for a better representation. 
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Table 2. Case D. Estimated position ^ (in units of Z, = 10^^ cm) 
and velocity v (in km s~') along the axis for the knots and output 
times of Fig.|9] 



apparent in the [S II] emissivity maps of Fig. |9] (left panels), at 



different output times of jet evolution. A secondary re-collimated 
jet appears at ^ = 130 L (see the frame at f = 460 yr), revealed 
by a light halo in the emissivity, corresponding to a high den- 
sity ring of matter. As already mentioned in Sect. 13.11 these re- 
collimation effects are not real, but arise from numerical effects 
due to the symmetry of the geometry. Five knots are visible to the 
left of the re-collimation point (knots A, B, C, D, E), and three 
to the right (knots F, G, H). Velocities can be estimated by look- 
ing at the temporal evolution of the emission function E{^), in 
the same figure (right panels). The position of each knot can be 
identified by tracking the spatial variations of the corresponding 
peaks in the 1-D plots. Results are shown in Table |2l The esti- 
mated velocities are seen to increase with the distance from the 
nozzle, from 3 km s"' to 64 km s"'. Even though these veloci- 
ties are still far from what observed in real jets (typically knots 
move at ^ 70% of the local flow speed), the trend of increasing 




Fig. 9. Case £). [S 11] emission map (left panels) and emissivity function E(^) (right panels) for different output times. 
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Radial velocity Vy versus R 
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Fig. 8. Case C. Radial velocity profile Vr(r) at t 
II L (solid line) and ^ = 21 L (dotted line). 
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velocities with distance from the source is invariably found in 
our simulations and it is a feature consistent with observations 
(lEisloffel & MundtlfT992h . 

A further comparison with observations of real jets can be at- 
tempted by measuring the decay in brightness of the knots over 
the beam length. If we measure this quantity for the first knots 
for the final output time t - 660 yr, their brightness decays with 
an estimated power law exponent a ^ -0.9. This value does 
not match observations and theoretical models, which foresee 
a X! -2. A va lue of g = -1. 9 has been found, as an exam- 
ple, in HH 30 (R aw et aLl figgel). In our simulation the value of 
a could be affected by the beam re-coUimation, which changes 
the brightness slope downward ^ ^ 130 L. This is why only the 
first four knots have been considered in the present estimate. We 
believe that such a problem will disappear in more realistic 3-D 
simulations. 



3.4. Restarting jets (case &) 

Case & retains the same parameters of the previous run. The dif- 
ference is that a temporal discontinuity has been applied on the 
initial conditions at the nozzle. Such time-varying initial con- 
ditions have nothing to do with the pulsating inflow conditions 
often invoked in the literature to generate moving knots. In case 
£ one makes the hypothesis that some discontinuity occurs at 
the jet source, namely that the flow at the nozzle is switched-off 
and that, after some time, it is switched-on again, so that for- 
mer initial conditions are restored. In this way a new, restart- 
ing jet forms, which bores its way in the wake of the old one, 
producing new, interesting features. The motivation for these in- 
flow conditions is provided by the observed formation of leading 
bow-shocks with temporal spacings of a few hundred years. In 
the present simulation a temporal variation has been set at the 
source of the jet with a time frequency small enough with re- 
spect to the intra-knot frequency, that one can still consider the 
overall structure as a steady jet. Note that this kind of variability 
is not meant to generate high-frequency internal working sur- 
faces as in the pulsating inflow models, since in our simulations 
knots form because of lOS along the jet beam. 

The results of the simulation for case & are displayed in 
Fig. [To] where [S II] 2-D maps and E{^) plots are shown at dif- 
ferent output times of the jet evolution. At f = 620 yr the source 
of the jet has been thus switched off. Both the inflow pressure 
and density, or mass injection rate, have been reduced by a fac- 
tor of 10. At f = 720 yr the original initial conditions have been 
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Table 3. Case £. Estimated position ^ (in units of L = lO''' cm) 
and velocity v (in km s ') along the axis for the knots and output 
times of Fig.fTOl 



restored, and a new jet forms at the source and starts to propa- 
gate in the wake of the old one. In the displayed plots we report 
emission maps (left panels) and E{£) functions (right panels), 
each row corresponding to different times of the evolution of 
this two-jet system, from t - 935 to f = 1050 yr This choice of 
output times allows us to study in detail the new jet, which ex- 
hibits a chain of knots quite c lose to the source, a typi cal feature 
observed in many YSOs (e.g.'Reipurth & Ballv" 200ll) . 

Interesting results are found by estimating the velocity of 
such knots (labeled A, B, C, D, E in Fig. [TOl ) from the po- 
sitions of the corresponding peaks in the integrated emission 
£'(^), as done in the previous sub-section. The kinematics of 
the new jet differs from the old one in two main aspects. The 
bow-shock propagates at a substantial fraction of the injection 
velocity Vjet - 200 km s"', while the external bow shock typ- 
ically propagates at a value which is about half of that speed. 
Moreover, the knots in the new jet are seen to move with rather 
high individual velocities, as reported in Table [3] This differ- 
ent behavior with respect to D is due to the fact that the new 
jet travels in the low density wake of the old one, rather than 
in the higher density ambient of the unperturbed ISM. The most 
remarkable results of this particular experiment are then the mor- 
phology and kinematics of the new jet knots, which appear close 
to the source and with relatively high-speed motions. 

4. Conclusions 

Our simulations show that under-expanded, light jets can nat- 
urally generate a pattern of emitting knots that possess proper 
motions, without invoking temporal variation of the source. In 
our scenario, knots are due to lOS, which are formed because 
of standard gas dynamical re-collimation processes, and their 
proper motion is due to the interaction with a highly time- 
dependent environment, namely the cocoon formed by the prop- 
agation of the jet head. The resulting knots are seen to survive 
radiative cooling and the synthetic images we obtained resemble 
qualitatively the observations of many HH objects: the individ- 
ual velocities are seen to increase with distance from the source 
and the knots' brightness, on the other hand, is found to decay 
over the beam length. However, the detailed properties of such 
knots, in terms of brightness, position in space, proper motion, 
and intra-knot spacing, obviously heavily depend on initial con- 
ditions. Exploring all, or even a large part, of parameter space is 
well beyond the scope and the possibility of this work, though 
some final considerations can be made: 

- the fact that steady inflow conditions can drive the formation 
of propagating emitting knots is a remarkable result in itself, 
since it has often been argued that steady jets could only form 
steady internal structures; 

- intra-knot spacing mostly depends on both 11, the pressure 
ratio and rjet, the nozzle radius. The correct choice of these 




Fig. 10. Case £. [S IIJ emission map and emissivity function E{^) at different output times. 
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parameters puts the knots length scale and the jet momentum 
loss in a realistic range; 

- the velocities of the knots can reach a significant fraction 
of the jet bow-shock propagation speed, which, in turn, 
is typically half of the velocity along the beam (basically 
the injection velocity). In restarting jets, simulating low- 
frequency (compared to the intra-knot frequency) variations 
of the inflow conditions, both the secondary bow-shock and 
the newly born knots are seen to propagate faster, reaching 
up to 40% of the local flow speed, which is not too far from 
what is observed in real jets (about 70%). This promising 
branch of numerical experiments has just been opened and 
will be more exhaustively explored in the future; 

- the need for 3-D calculations arises from some significant 
discrepancies between numerical results and obervationally 
determined properties, such as the emissivity decaying ex- 
ponent with the distance from the source. We claim here that 
such differences could arise from artificial re-collimation ef- 
fects due to the assumed hypothesis of axisymmetry, but 3-D 
simulations are needed to prove this statement. 

In conclusion, our results show that lOS provide a natu- 
ral, efficient mechanism for the formation of radiatively emit- 
ting knots which possess most of the observed features, such 
as proper motions with increasing velocities along the jet beam. 
However, claiming that they are the only driving mechanism is 
not realistic, and most probably lOS work in co-operation with 
other mechanisms. In this framework, knots arising from lOS 
or from local working surfaces generated by inflow conditions 
fluctuations could either co-exist or work separately in different 
objects. 
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